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A flow of electrically conducting fluid in the presence of a steady magnetic field has 
a tendency to become quasi two-dimensional, i.e. uniform in the direction of the mag- 
netic field, except in thin so-called Hartmann boundary layers. The condition for this 
tendency is that of a strong magnetic field, corresponding to large values of the dimen- 
sionless Hartmann number (Ha >> 1). This is analogous to the case of low Ekman 
number rotating flows, with Ekman layers replacing Hartmann layers. This has been at 
the origin of the homogeneous model for flows in a rotating frame of reference, with 
its rich structure: geostrophic contours and shear layers of Stewartson [Ste57, Ste66], 
Munk [Mun50] and Stommel [Sto48]. In magnetohydrodynamics, the characteristic sur- 
faces introduced by Kulikovskii [Kul68] play a role similar to the role of the geostrophic 
contours. However, a general theory for quasi two-dimensional magnetohydrodynamics 
is lacking. In this paper, a model is proposed which provides a general framework for 
quasi two-dimensional magnetohydrodynamic flows. Not only can this model account for 
otherwise disconnected past results, but it is also used to predict a new type of shear 
layer, of typical thickness Ha -1 / 4 . Two practical cases are then considered: the classical 
problem of a fringing transverse magnetic field across a circular pipe flow, treated by 
Holroyd and Walker [HW78] , and the problem of a rectangular cross-section duct flow in 
a slowly varying transverse magnetic field. For the first problem, the existence of thick 
shear layers of dimensionless thickness of order of magnitude Ha -1 / 4 explains why the 
flow expected at large Hartmann number was not observed in experiments. The second 
problem exemplifies a situation where an analytical solution had been obtained in the 
past [WL72] for the so-called "M-shape" velocity profile, which is here understood as an 
aspect of general quasi two-dimensional magnetohydrodynamics. 



1. Introduction 

The relative influence of Lorentz forces versus viscous forces is measured by the dimen- 
sionless Hartmann number Ha. At large Hartmann number, the flow can be analyzed in 
terms of an inviscid core flow and viscous shear layers. This is very similar to the case 
of low Ekman number flows in a rotating frame of reference, for which two-dimensional 
models are commonly used. The tendency of flows of electrically conducting fluids under 
a DC magnetic field to become two-dimensional has been known since the work of Hart- 
mann and Lazarus [HL37a, HL37b]. A number of analytical [She53, CL61, Hun65] and 
experimental [Mur53, Bra67] studies have confirmed this trend in the laminar regime and 
even in the next stage of so-called two-dimensional turbulence [Leh55, Kol72]. It was then 
natural to take advantage of the well-known structure of the Hartmann layers to build 
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a model for the two-dimensional flow. This was done by Sommeria and Moreau [SM82], 
in terms of the core velocity and pressure. There have been also other two-dimensional 
models using the electrical potential and pressure as main variables [MB94]. However 
these models have all in common that the flow is supposed to be contained between 
two parallel planes with a uniform transverse magnetic field. This constitutes a very 
particular case as the subsequent developments show. 

The work of Kulikovskii [Kul68] brought more generality since he considered electrically 
insulating cavities of arbitrary shapes in arbitrary non-uniform magnetic fields. He shows 
that there is still a kind of two-dimensionality as the variation of physical quantities like 
pressure or electric potential along the magnetic lines can be obtained from the equations 
of the core of the flow. In addition, he shows that the quantity J ds/\\B\\, integral on a 
magnetic line of the inverse of the magnetic field intensity, plays a crucial role. Surfaces 
made of magnetic lines with the same value for this integral are called characteristic 
surfaces: Kulikovskii shows that the electric current and velocity have a tendency to 
lie on these characteristic surfaces at large Hartmann number. These surfaces can be 
considered as surfaces of least resistance to the flow. A flow can be forced to cross them 
but this will cost more in terms of pressure gradient required than in the case of a flow 
along themf. In a second paper [Kul73], Kulikovskii applies his asymptotic structure to 
particular configurations, showing the powerful predictions that can be made using the 
concept of characteristic surfaces. 

Later, the ideas of Kulikovskii have been applied by Holroyd and Walker [HW78] in the 
case of a fringing transverse magnetic field to a circular duct. The transverse magnetic 
field is uniform on both sides of a step change. The authors follow Kulikovskii and assume 
that in the non-uniform region, the flow follows strictly the characteristic surfaces and 
that on each side it would then come back towards a fully-developed regime on a much 
larger length scale Ha 1 / 2 , as shown from scaling analysis!- They then re-derive a two- 
dimensional model with uniform magnetic field but non-uniform depth of the cavity 
in the direction of the magnetic field, in terms of pressure and electric potential. This 
model is then solved using expansions in eigenfunctions on both sides, where pressure 
and potential are matched in the central non-uniform magnetic field region. In doing so, 
they obtain the large Hartmann number asymptotic limit of the duct flow with fringing 
magnetic field. Unfortunately, experiments performed under a large Hartmann number 
by Holroyd [Hol79], Ha = 523 based on the radius of the pipe, showed that the observed 
flow was significantly different from the predicted one: it is modified in the same way 
as expected but much less spectacularly The reason for this discrepancy could not be 
explained satisfactorily by invoking possible viscous effects in the core or inertial effects 
since the Hartmann and interaction numbers were both large enough to discard such 
effects. 

In a second attempt, Hua and Walker [HW89] derive two-dimensional equations not 
only for variable depth but also for a non-uniform magnetic field. The assumption of 
so-called parallel magnetic lines is made, whereby the magnetic field is supposed to be 
pointing nearly in the same direction (z), although its intensity depends on the other 
two directions (x and y). This model could be solved numerically without the need to 
match the pressure and potential in the non-uniform magnetic field region. The results 
were consistent with experimental observations at comparable values of the Hartmann 

f For the same velocity, the ratio of pressure gradient required to drive a flow across charac- 
teristic surfaces is of the order of a factor Hartmann times the pressure gradient necessary to 
drive a flow along these surfaces. 

J This scaling appears for the first time in [WL74]. 
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number. Moreover, the numerical solution could be obtained for higher values of the 
Hartmann number up to 7000. The authors remark that, even at this large value of 
Hartmann number, the computed flow is far from the asymptotic prediction of Holroyd 
and Walker. As they expect the asymptotic regime to be reached when the square root 
of the Hartmann number is large compared to unity they conclude rather unconvincingly 
that V7000 ~ 84 is not large enough compared to unity to approach the asymptotic 
regime. 

A closely related problem is that of a cylindrical duct with expansion in a uniform 
transverse magnetic field [WL74] . The topology of the characteristic surfaces is identical 
to that of the duct of constant diameter with varying transverse magnetic field. This 
study is purely analytical, based on matching solutions on both sides of the region of 
varying diameter. 

Another important case where the characteristic surfaces play a crucial role is that of 
the so-called entrance problem for rectangular duct flows. When the flow enters a region 
of varying transverse magnetic field, the flow appears to concentrate along the sides of 
the duct parallel to the magnetic field, forming a so-called M-shape velocity profile. This 
was reported from experimental observations in a range of papers [SSS70, BS66, Hol79, 
Hol80]. For this configuration again, Walker and co-workers played a major role in the 
analysis of the flow [WLH72, WL72, SW99]. 

Throughout this paper, the question of the relevance of a two-dimensional MHD model 
will be addressed. It is hoped that the reader will be convinced that two-dimensional 
equations of the type of those presented in [HW78, HW89] do represent MHD flows 
correctly, but that their scaling analysis has not received enough attention. Such analysis 
provides a global understanding of MHD flows. In particular, it will be shown how some 
known results can be rederived (e.g. the M-shapc profile) and how an open question 
arising in the case of a fringing transverse magnetic field to a circular duct can be clarified. 
Also, the comparison with the case of rotating flows will be made. As the two-dimensional 
analysis is more advanced and better understood for rotating flows, this comparison will 
hopefully give support to the proposed MHD analysis. In addition, strong magnetic fields 
are becoming widely available for research and industrial purposes, and understanding 
the asymptotic limit of large Hartmann number flows is becoming increasingly useful 
in metallurgy or other processes involving liquid metals, like liquid metal cooling or 
spallation neutron sources. 

2. Three-dimensional asymptotics 

In this section the different types of three-dimensional boundary layers appearing in 
MHD flows or flows in rotating systems of reference are reviewed briefly. These well-known 
features are recalled here for future comparison in section 3 with typical two-dimensional 
structures and to stress the analogy between MHD and rotating flows: the theory of 
rotating flows is more advanced than that of MHD flows and the main point of this 
paper is to show that the diversity of scaling laws found in the homogeneous model of 
rotating flows also exists for MHD flows. 

2.1. MHD three-dimensional asymptotics 

The dimensionless vector position x, magnetic field B, velocity field u, electric current 
density j, pressure field p and electric potential field <j) ar( 3 derived from their corre- 
sponding dimensional quantities using the scales H, B , v/H, voB^jH, pv 2 / H 2 and 
vBq respectively, where H is a length-scale of the configuration, Bq a typical value of 
the applied magnetic field, p the density of the fluid, v its kinematic viscosity and a 
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its electrical conductivity. The magnetic field is imposed externally and the magnetic 
Reynolds number is assumed to be small so that the dimensionless magnetic field B is 
divergence- free and curl- free according to Maxwell's equations. Moreover, inertia is also 
neglected so that the dimensionless expressions of Navier-Stokes equation and Ohm's law 
take the following form: 

= -Vj) + ffa 2 jxB + V 2 u, (2.1) 
j = -V^ + u x B. (2.2) 

In addition, the velocity field u and electric current density j must be divergence-free 
and on the boundary of the cavity enclosing the fluid, the velocity must be zero as well 
as the component of the electric current density parallel to the normal vector, since only 
electrically insulated cavities are considered in this paper. One could consider a field 
of volume forces in the Navier-Stokes equation, like buoyancy forces in the Boussinesq 
approximation, without affecting the analysis of the structure of the flow, provided this 
force field does not vary on a short length-scale. In the set of equations (2.1) and (2.2), 
the Hartmann number Ha = y/ 'a j pvB^H is the single dimensionless parameter. 

Asymptotic analysis refers to the limit of infinite Hartmann number Ha. Since Hart- 
mann, it has been known that the flow can be divided into so-called core regions of finite 
size and layers of small thickness. A simple way to put these structures in evidence is to 
perform a local analysis of the governing equations (2.1), (2.2), so that one can consider 
that the magnetic field is uniform. Let us denote z the local coordinate aligned with 
the magnetic field and x, y other coordinates such that (x, y, z) is a direct orthonormal 
system of reference. Taking twice the curl of (2.1) and substituting V x j using the curl 
of (2.2) leads to the following equation for the velocity field: 

H ^q^-(V 2 ) u = 0, (2.3) 

where Hai = Ha ||B|| represents the local Hartmann number. The electric current density 
j can be shown to obey the same equation. The three-dimensional structures of MHD 
flows can be derived from the differential operatorf Ha 2 d 2 /dz 2 — (V 2 ) . 

2.1.1. Core regions 

When the Hartmann number tends towards infinity, the bi-Laplacian term cannot 
compete in finite-size regions and one gets a region where the second derivative of the 
velocity vanishes. In the absence of a strong electric current along the magnetic lines, it 
can be shown that a two-dimensional velocity field develops in these regions of finite ex- 
tent. Even in the general case of a non- uniform magnetic field, a kind of two-dimensional 
state exists, since the variations of velocity along the magnetic lines are readily obtained 
from the governing equations. 

2.1.2. Hartmann layers 

Where the magnetic field is not parallel to the wall, a Hartmann boundary layer 
develops. Obviously, the direction of maximum variation is the normal direction to the 

f For completeness, it is recorded that this operator can be decomposed into two operators as 
follows: Haf d 2 /dz 2 — (V 2 ) = {Hai d/dz — V 2 ) o {Hai d/dz + V 2 ) , where each of these opera- 
tors governs its corresponding Shercliff variable [She53] . These variables are linear combinations 
of the velocity and (induced) magnetic field, although they are a different linear combination 
than Elsasser variables [Els50] . One must introduce another dimensionless number, the magnetic 
Prandtl number Prm = fiau, where /i is the magnetic permeability. While Elsasser variables 
are u ± A, with A the Alfven velocity B/^/p/I, Shercliff variables are u ± Prm -1 / 2 A. 
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wall n and equation (2.3) can be simplified into: 

where n is the distance to the wall. Hartmann layers are solutions to this equation and 
the tangent velocity varies exponentially in them on a typical length-scale Ha -1 (B.n) -1 . 

2.1.3. Parallel layers 

Viscous terms represented by the bi-Laplacian term in (2.3) can also play a role in 
thin regions (boundary or free-shear layers) containing the magnetic field direction. In 
this case, equation (2.3) takes the following form: 

where n is again a coordinate perpendicular to the layer. The solutions to this equation 
are parallel layers. It is difficult to obtain a general solution in these layers, but a scaling 
analysis shows that if these layers stretch along a distance of order unity along the 

— 1/2 

magnetic field lines, their thickness must be of order Ha l . Parallel layers were first 
analyzed by Shcrcliff [Shc53]. 

Interestingly, another result can be derived from equation (2.5) in the case when the 
cavity is of infinite length in the direction of the magnetic field. This gives a limit to the 
size of the core region where the flow is two-dimensional. There can be variations along 
the magnetic lines on a large length-scale: if the typical cavity length scale perpendicular 
to the magnetic lines is of order 1, both terms in (2.5) can balance provided the length- 
scale along the magnetic lines is of order Hai. This situation arises typically when an 
object of size of order unity is placed in an infinite region: its effect on a flow around it 
stretches along the magnetic lines as far as a distance of order Ha. 

In the case of parallel layers and when the magnetic field is non uniform, the local 
analysis is not rigorous as one has to consider variations of the velocity on a length- 
scale of order unity or more along magnetic lines. Equation (2.5) is therefore not correct. 
However the scaling analysis remains valid: in a paper by Todd [Tod68], the rigorous 
analysis of parallel layers developing along curved magnetic lines is performed and it is 
shown that these layers are still lying on magnetic lines and arc still of typical thickness 
Ha- 1 ' 2 . 

2.1 A. Roberts layers 

Roberts layers develop as a curved wall approaches a point where it becomes tangent 
to the magnetic lines. The Hartmann layer solution presents a singularity as its thickness 
diverges. The divergence is resolved by a so-called Roberts layer [Rob67]. If the curvature 
is of order unity at the point of contact, a small region elongated in the magnetic field 
direction must satisfy the geometric constraint that its thickness across magnetic lines 
scales as the square of its length along these lines. By scaling analysis of equation (2.5) 

2/3 1/3 

the only possibility is a thickness of order Ha l and a length of order Ha t along 
the magnetic lines. 

2.2. Rotating flows three-dimensional asymptotics 

Now, instead of applying a steady magnetic field B, it is supposed that the reference 
system considered is rotating with the rate Q along the z axis. A Coriolis force replaces 
the Lorcntz force in the momentum equation. The same dimcnsionlcss variables are used 
for distance, velocity and pressure as in section 2.1. The governing equation can then be 
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Figure 1. Sketch of three-dimensional structures in MHD and rotating flows: (a) Hart- 
mann / Ekman layers and MHD parallel layer / inner Stewartson layer; (b) Roberts layer / sin- 
gular Ekman layer; (c) MHD wake / 'Taylor' column. 

written as follows: 

= -Vp + 2£r 1 u x e z + V 2 u, (2.6) 

where E = i//(QH 2 ) is the dimensionless Ekman number. Taking the curl of equation 
(2.6) leads to an equation involving u alone: 

= 2E- 1 |f + V 2 (V x u). (2.7) 

One can even differentiate this equation with respect to z to eliminate the curl operator 
and obtain a form suitable for analytical investigations: 

= 4iT 2 §^ + (V 2 ) 3 u. (2.8) 

Under the form (2.7) or (2.8), one can extract the possible structures for the velocity 
field using scaling analysis (see [Gre68] for a classical treatment). When a boundary layer 
develops on a wall that does not contain the direction of rotation, a layer of dimensionless 
thickness E 1 / 2 develops: this is an Ekman layer and corresponds to the Hartmann layer 
in MHD. It can also be seen that shear layers containing the direction of rotation can 
develop. If the dimensionless length of the cavity in the direction of the rotation is of 
order one, these layers must have a thickness of order E 1 ^ 3 : they are the strict equivalent 
to the MHD parallel layers. When an object of dimensionless size one is held in place 
in an infinite domain, then its 'wake' extends a distance E~ x in the direction of the 
rotation. Finally, similar regions to Roberts layers develop at the point of a curved 
wall where it becomes tangent to the direction of rotation: these regions have a typical 
length-scale of order E 1 / 5 and E 2 / 5 in the direction and perpendicular to the direction of 
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rotation respectively, if the curvature is of order unity. They correspond to a singularity 
of diverging Ekman layers. 

At this point, it can be seen that there is a close correspondence between the possible 
structures in MHD and in rotating flows according to the three-dimensional equations 
(see Fig. 1 for a pictorial summary). However, it is not clear how one can possibly extract 
the typical E 1 / 4 main component of Stewartson's layers from the scaling of the differential 
operator AE~ 2 d 2 /dz 2 + (V 2 ) 3 appearing in equation (2.8). This becomes obvious when 
two-dimensional equations are derived. 

3. Two-dimensional models 

A cavity of dimensionless length-scale unity in the direction of the imposed magnetic 
field or rotation is considered. As we have seen above in section 2, the largest three- 
dimensional structures in the direction perpendicular to the magnetic field or rotation 
direction are of order Ha" 1 ! 2 or E 1 / 3 . Hence, if only larger perpendicular scales are 
considered, the flow consists simply of a two-dimensional core bounded by Hartmann or 
Ekman layers. This is at the origin of the two-dimensional models. 

3.1. MHD two-dimensional model 

In most cases, MHD core flows are two-dimensional, except when the forcing and bound- 
ary conditions are such that an odd velocity profile is generated along a magnetic field 
line (a linear variation in the core). In this case, the MHD braking is very strong: this 
is the case referred to as 'singular symmetry' in [AGM96]. The other symmetry called 
'regular symmetry' corresponds to an even velocity profile (which is nearly uniform in 
the core) along magnetic lines, and the MHD braking is Ha times weaker. Hence, the 
general case bears resemblance with the 'regular symmetry' and corresponds indeed to 
the idea of a two-dimensional flow. 

In the following analysis a cavity of general shape and a non-uniform magnetic field 
are considered. Following most authors, the approximation of straight magnetic lines 
is made here. In this approximation, the curvature of the magnetic lines is ignored, 
the intensity of the magnetic field is taken as constant on each magnetic line but can 
vary from one line to another. This is clearly an approximation intended to simplify 
the analytical calculations, but it has been shown to be safe in the few studies where 
this approximation was not made [Kul68, Tod68, AGM96]. This approximation does not 
affect the qualitative features of the flows, although the imposed magnetic fields are not 
physical since they are rotational. 

As shown in Fig. 2, the geometry of the electrically insulated cavity consists of the 
space situated between two surfaces, an upper surface S u and a lower surface S l , defined 
in a orthonormal coordinate system (x,y,z) by the two functions z u (x,y) and z l (x,y) 
respectively. The z axis is chosen to coincide with the magnetic field direction, in the 
straight magnetic lines approximation. Its intensity depends on x and y, B = B z (x, y)e z . 
The functions z u , z l and B z arc dimensionless functions of the dimensionless coordinates 
x and y, where the arbitrary scales H and Bq for length and magnetic field intensity 
defined in section 2 continue to be used. 

It is assumed that, in the core of the flow, the velocity and electric current density 
components perpendicular to the magnetic field are independent of z. Using for instance 
the value of these fields on a reference plane z = 0, one can define the two-dimensional 
fields: 

u = u 0x (x,y) 
_ u 0y {x,y) 



Jo = 



3ox(x,y) 
hy(x,y) 



(3.1) 
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Figure 2. Geometry and magnetic field for two-dimensional analysis 



Pressure and electrical potential are also independent of the magnetic field direction z 
in the core, and their two-dimensional representatives are denoted po and <f>o. Consider- 
ing the restricted versions of the general three-dimensional equations (2.1) and (2.2) in 
the plane (x,y) in the core, two-dimensional equations governing Uo, jo, Po and 4>q are 
obtained: 

O = -X7p o + Ha 2 B z j xe 2 + V 2 u , (3.2) 
jo = -V</>o + B z u x e z . (3.3) 

These governing equations are not sufficient to describe the MHD two-dimensional flows. 
First, one needs to take into account the global conservation of mass and of electric 
charge. Second, one needs to take into account the contribution due to the Hartmann 
layers in these conservation relationships. 

The mass and charge conservation laws will be expressed in terms of global two- 
dimensional mass and electric charge flux densities Q and I, taking into account the 
core and Hartmann layers contributions. The core contribution is obtained by multi- 
plying the two-dimensional vector fields uo and jo by the length of the magnetic line 
within the cavity z u — z . The Hartmann contribution to the mass flux is neglected as 
this corresponds simply to a deficit, as measured by the displacement thickness of the 
Hartmann boundary layer. It introduces an error of order i?a -1 , comparable to the ap- 
proximations already made when assuming a two-dimensional core flow. On the contrary, 
the Hartmann layer contribution to the electric charge flux is fundamental in MHD. It 
is proportional to the core velocity adjacent to the Hartmann layer and in the direction 
perpendicular to both this core velocity and the wall normal vector. The core velocity 
adjacent to the Hartmann layer u c must be tangent to the wall and the integral of the 
electric current over the Hartmann layer thickness I# a depends on the wall normal unit 
vector n and u c (sec [HS71]): 

Ijj = — sign(B.n) Ha^ 1 u c x n. (3.4) 

This surface electric current density needs to be expressed in terms of the coordinates x 
and y which will be used in the following two-dimensional analysis. To this end, let us 
calculate the amount of electric current flowing across a small line clement spanned by 
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Figure 3. Local tangent wall surface and relationship between the core velocity and total 
electric current flowing in the Hartmann layer 



dx and dy (see Fig. 3). This line element corresponds to a line element dl on the tangent 
plane to the surface and, from the components uq x and u oy of uo, one can also build the 
tangent core velocity: 



dl 



dx 
dy 

dx^f- 

OX 



dyl 



UOx 

UOx-Qx- 



u oy-dy- 



(3.5) 



where the upper surface S u is considered here, the treatment for the lower surface being 
similar. The amount of electric current flowing across the line element dl is proportional 
to the magnitude of the vector dl x (u c x n). This vector can be shown to be equal to 
(c£l.u c ) n. Hence, with a factor Ha -1 , the flux of electric current through dl is given by: 



dx 



UQx + u 0x 



dz u 
dx 



Uoy 



dz u dz w 
dx dy 



dy 



UOy + U y 



dz u 
dy 



UQx' 



dz u dz u 
dx dy 



(3.6) 



The two-dimensional vector flux T^ a in the (x, y) coordinate system that represents that 
electric current flux has components: 



TU 

1 Hax 

TU 

1 Hav 



= Ha 



-U y ~ U y 

i n . /(Jz"^ i dz u dz^ 

U 0x +U 0x {- 5F ) +U y- d - x -- d - y - 



dz u gz" 
v x dx dy 



(3.7) 



The lower surface produces a similar electric current flux. Both contributions are ex- 
pressed by a two-dimensional tensor, T u and T l for the upper and lower surfaces respec- 
tively: 



l Ha 



with 



T" 



dz u dz u 
dx dy 



dz u dz u 
dx dy 



_dz^d£_ _j 
dx dy 

V ox J ox oy 



(3.8) 



(3-9) 



It is now possible to express the total mass and charge two-dimensional fluxes, including 
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the core and Hartmann layers contributions: 

Q = (z u -z l )u , I=(z u -z% + I u Ha + I l Ha . (3.10) 

For impermeable and electrically insulated surfaces, these two-dimensional flux densities 
arc divergcnccless, hence can be expressed in terms of the streamfunctions ip and h. 

Q = Vi/' x e 2 , I = Whxe z . (3.11) 

Taking the curl of the two-dimensional equations (3.2) and (3.3) eliminates pressure and 
electric potential fields and yields the following equations, which have only one component 
along z: 

= -Ha 2 [V.(BJo)] + V 2 (V x u ) , (3.12) 
Vxj = -[V.(B 2 u )]e z . (3.13) 

These equations can then be expressed entirely in terms of ip and h, using (3.11), (3.10) 
and (3.8) for substitution. For the sake of concise expressions, let us introduce d = z u — z l 
the depth of the cavity in the direction of the magnetic field and T = T u + T l the sum 
of the upper and lower Hartmann electric currents in terms of the core two-dimensional 
velocity. During the substitution, the contribution of these Hartmann currents to equation 
(3.13) is neglected, while their contribution to (3.12) is retainedf: 



(3.14) 



t)-'(t)"* 

Both equations have almost a symmetrical structure. The curl of Navier-Stokes (3.14) 
has a diffusion term for ip on the left-hand side: it is related to the two-dimensional 
vorticity of the flow. The equation has a term of the form VK x Wh on the right- 
hand side, with K = B z /d. This fundamental term corresponds to the flow of total 
electric current crossing the characteristic surfaces. Finally, the last bi-Laplacian term 
on ip represents the effect of viscosity in the two-dimensional flow. This last term will 
be shown to be negligible in most cases, and then equation (3.14) expresses that when 
the electrical current goes across characteristic surfaces, this creates vorticity. Equation 
(3.15) is similar, except that there is no equivalent to the last term of (3.14). It says 
that when the flow goes across characteristic surfaces, this creates a region of curl for the 
electric current. 

It is important to stress that the characteristic function K = B z /d appears naturally 
in these two-dimensional equations, as could be expected from previous analyses [Kul68, 
AGM96]. The derivation and results presented above are close to those presented in 
[AlbOl] . This analysis is also very close to the analysis by Holroyd and Walker [HW78] and 
Hua and Walker [HW89] , where the main variables are the pressure and electric potential 
instead of the streamfunctions ip and h considered here. Also, they have neglected from 
the beginning the viscous friction in the two-dimensional flow. Although the analysis 
above about the significance of the terms is not carried out in papers [HW78, HW89], it 
can be checked after rearranging different terms that equations (17) and (18) in [HW89] 
display the same structure as (3.14) and (3.15) of this paper, with the characteristic 
function B z /d playing an important role. 

Equations (3.14) and (3.15) are going to be solved rigorously in sections 4 and 5. 

t The motivation for these different treatments is that the Hartmann currents are nearly 
curl- free while their divergence can be significant. 
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Nevertheless, it is important to perform a scaling analysis of these equations to identify 
the possible structures that can be expected in two-dimensional flows. For the purpose of 
the scaling analysis, the depth d and magnetic field B z can be considered to be uniform 
and unity in all terms except for the two terms involving the gradient of K , since variation 
of K = B z /d is crucial. Moreover, in the spirit of a local analysis, a change of coordinates 
is introduced, in which (x,y) is replaced by orthonormal intrinsic coordinates (r, s), 
where r varies in the direction perpendicular to the characteristic surfaces, while s varies 
along the characteristic surfaces. The symbol G is chosen to denote the magnitude of 
the gradient of K: G = \\V(B z /d)\\. With these simplifying assumptions and notations, 
equations (3.14) and (3.15) take the form: 

Bh 

HaW 2 tp = Ha 2 G— + (W 2 ) 2 ^, (3.16) 
os 

V 2 h = G^ (3.17) 

os 

Taking the Laplacian of (3.16) and substituting h using (3.17) provides the fundamental 
differential equation governing the two-dimensional streamfunction: 

ifc(V 2 ) 2 ^ifc 2 G 2 ^$ + (V 2 )V (3.18) 
os z 

Under this form, it is convenient to perform a scaling analysis of the two-dimensional 
MHD flows. First of all, in the absence of thin layers, the term with the largest power of 
Ha is dominant. Consequently the second variations of ip along the characteristic surfaces 
must be zero. In fact, very often the first variation of ip will be zero along theses surfaces, 
in order to minimize Ohmic dissipation (see equation (3.17)). This corresponds to the 
'core regions' of two-dimensional flows. Different type of layers can appear, depending 
on whether they develop in a direction parallel to characteristic surfaces or not. If they 
develop along a wall not parallel to the characteristic surfaces, then equation (3.18) 
becomes: 

Hap^=Ha 2 G 2 co S 2 0p^ + p^, (3.19) 
cm 4 on z on° 

where n is the normal direction to the wall and 9 is the angle between the normal to 
the wall and the direction of characteristic surfaces. The thinnest structure that can 
emerge from equation (3.19) is obtained when the last term (viscous core friction) is 
balanced with the first term (Hartmann layer friction). This leads to a typical thickness 
of Ha -1 / 2 and must immediately be discarded as this scale corresponds to the scale of 
the three-dimensional parallel layers. Indeed, the three-dimensional analysis in section 
2.1 shows that for a length scale of order Ha~ x l 2 or less in the direction perpendicular to 
the magnetic field, the variations of velocity in the core of the flow arc of order unity in 
the direction of the magnetic field: hence the two-dimensional model ceases to be valid. 
This means that the core friction will always be negligible compared to Hartmann layer 
friction at larger scales. The other possibility is to balance Hartmann friction with the 
topographic term: this leads to a layer of typical thickness Ha~ 1 / 2 G~ 1 cos~ 1 6. Again, if 
G and cos# are of order unity, this two-dimensional layer is not physical as it is replaced 
by a Ha~ x l 2 three-dimensional parallel layer. However, if GcosO is very small compared 
to unity, this thickness Ha~ 1 / 2 G~ 1 cos~ 1 6 is larger than the parallel layers and consti- 
tutes a valid two-dimensional layer. This type of layer has been analyzed by Walker and 
Ludford [WL74] who started from the three-dimensional equations. These layers will be 
examined in section 5 from the easier point of view of the two-dimensional equations. 
They are similar to Stommel layers [Sto48] in rotating flows, as the main balance involves 
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topographic effects and friction in Ekman layers. However, as will be seen in the next 
section, the Stommel solution does not apply, as the core viscous friction turns out to be 
larger than the Ekman friction in rotating flows. Hence, Stommel layers arc replaced by 
thicker Munk layers [Mun50]. 

If one considers now a layer developing along a characteristic surface, equation (3.18) 
can be written: 

Again, it can be seen that the last core friction term will always be negligible on scales 
larger than Ha -1 / 2 and can thus be ignored in this two-dimensional analysis. 

The scaling analysis depends on the scale along the characteristic surface. If this scale 
is of order unity (d/ds ~ 1), then the balance between the Hartmann friction term and 
the topographic term leads to a layer of typical thickness Ha~ x ^G~ x l 2 which is much 
thicker than parallel layers even if G is of order unity. This type of layer will be discussed 
in section 4. 

Finally, if the cavity is infinite in the direction of the characteristic surfaces and of 
size unity in the perpendicular direction, equation (3.20) can provide the typical long 
length-scale of development of structures along characteristic surfaces. If d/dr ~ 1, 
then the balance between Hartmann friction and topographic effects leads to a scale 
Ha x / 2 G. Walker and co-workers [WL74] had found that a distance of order Ha 1 / 2 was 
the length-scale necessary from the position where the magnetic field was varying to 
reach a fully-developed flow in a pipe with transverse magnetic field. 

There can be other possibilities like the case when a wall would become tangent to the 
characteristic surfaces: similarly to Roberts layers, they would correspond to a singularity 
in the Ha~ 1 / 2 G~ 1 layer. For a curvature of order unity, scaling analysis of (3.20) when 
the length-scale perpendicular to the wall scales as the square of the length-scale along 
the characteristic surfaces, one gets a region of size HaT^G -2 ' 3 by Ha~ 1 l % G~ 1 / 3 . The 
structures discussed above are believed to be the most typical ones and are illustrated in 
Fig. 4. 

As discussed above, the core viscous term (V 2 ) ip in equation (3.16) is always negli- 
gible when the two-dimensional equations are valid. Without this term, equation (3.16) 
may be multiplied by Ha~ 3 ^ 2 , then successively added to and subtracted from equation 
(3.17) providing the following equations: 

V 2 (i>/Ha 1 ' 2 ±h)= ±Ha l ' 2 G^- s (i'/Ha 1 / 2 ± h) . (3.21) 

Those familiar with Shercliff variables will notice that ip/Ha 1 ^ 2 ± h play here a corre- 
sponding role. The privileged direction is that of characteristic surfaces instead of the 
direction of the magnetic field for Shercliff variables. This may help to accept the typical 
length-scales derived above. This also indicates that the magnitude of the electric current 
(h) is of order HaT 1 ! 2 that of the velocity (ip) in general. 

3.2. The homogeneous model of rotating flows 

The corresponding two-dimensional model in the case of rotation is called the the homo- 
geneous model. Its derivation can be found in many textbooks (e.g. Greenspan [Gre68]) 
and is sketched here in order to emphasize the similarity with the MHD model. At the 
upper and lower surfaces, the role of the Ekman layer is now to produce a cross-flow with 
respect to the core flow direction. The integral of this cross flow scales as E 1 ! 2 times 
the core velocity. For this term and for the purpose of this approximate derivation, the 
angles between the direction of rotation and the wall normals are neglected. The total 




Figure 4. Two-dimensional MHD structures. 

vertically-integrated volume flow rate can then be expressed: 

Q = du Q + E 1/2 e z x u . (3.22) 

The core two-dimensional equation is now obtained by restricting equation (2.6) to the 
(x, y) components in the core of the flow: 

= -Vp + 2£ , ~ 1 u x e^+V 2 u . (3.23) 

Its curl takes the form: 

0=-2r 1 (V.u )e 2 + V 2 (Vxu ). (3.24) 

Again, one can use the streamfunction tp, defined in (3.11) and substitute it for Uo in 
(3.24). The Ekman layer cross-flow is neglected in the viscous term and plays a crucial 
role in the other term. This results in the homogeneous model: 

2 E-W V. (?f) e z =2^V Q) x V V , + V 2 (v. (?f) ) e z . (3.25) 

A similar change of two-dimensional coordinates to that in the previous section is 
operated, where the new (r, s) system is defined such that r is constant along lines of 
constant depth, the gcostrophic contours, and s is changing along these lines. Defining 
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(3 = ||V(l/<i)|| in accordance to the so-called /3-plane approximation - !", equation (3.25) 
can be written locally: 



under the assumption that d is of order unity. Although the order of derivation of the 
terms and the powers of Ekman number are different from those of Hartmann number 
in the MHD case, there is a strong analogy with the MHD equation (3.18). The term on 
the left-hand side represents the contribution of the Ekman (rcsp. Hartmann) layer, the 
first term on the right-hand side corresponds to the effect of topography while the last 
term is related to shear forces within the two-dimensional bulk flow. 

In terms of boundary-layer structures, equation (3.26) can be used to estimate the 
thickness of possible boundary layers arising in rotating flows. For a cavity of constant 
depth, there can be an equilibrium between the first and last term in (3.26) on a typ- 
ical length-scale E 1 / 4 . This is the thickness of shear boundary layers between regions 
of different core velocity [Ste57]. Along western boundaries, the equilibrium between 
the second and third terms brings a length-scale E 1 / 3 / '/3, corresponding to Munk layers 
[Mun50, CW77, Wal75]. The equations are not symmetrical between East and West, 
owing to the term fid^p/ds in equation (3.26), and no boundary layer can develop on 
eastern boundaries. The development length of a flow along geostrophic contours after a 
disturbance can be estimated as E' 1 / 2 [Wal74]. 

4. Circular duct with fringing magnetic field 

A long electrically insulating circular duct is submitted to a transverse magnetic field 
(Fig. 5). The axis of the pipe is taken to be the x-direction, while the direction of 
the transverse magnetic field is taken to be the z direction. All dimensions are made 
dimcnsionless using the diameter of the pipe. From x = — oo to x = —2 the magnetic 
field intensity is assumed to be uniform. Its value Bq is used as our scale for magnetic 
field, hence its dimcnsionless value is unity, while from x = 2 to x = oo, the dimcnsionless 
magnetic field is uniform with an intensity twice as large. Between x = —2 and x = 2, 
the dimensionless magnetic field is assumed to have the form 



which makes B z continuous, as well as its first and second derivatives, from x = — oo to 
x = oo. According to the straight magnetic field lines approximation, the other compo- 
nents of the magnetic field are assumed to be zero, although this field does not satisfy 
the condition of a curl-free magnetic field strictly speaking. 

The characteristic surfaces contain the z direction and are best made visible in the 
{x,y) plane where they appear as lines. On Fig. 6, some characteristic surfaces are drawn: 
they correspond to the case of a circular cylinder and to the distribution defined by the 
magnetic field (4.1). However, the x and y coordinates are stretched differently. The 
figure is drawn for the particular case of a duct of length 8, but this is irrelevant since 
the characteristic surface are independent of x when x < — 2 or x > 2. In fact, the 
length of the cavity used for numerical calculations will be adapted to each value of the 
Hartmann number (see section 4.1) so that it is always longer than the development 
length of the flow. 

f The traditional /3-plane approximation corresponds to the linearized variation of latitude, 
but this is equivalent to a change in depth, as far as the two-dimensional model is concerned. 



2 E- 1 ' 2 VV - 2 E" 1 (3^- + (V 2 ) V, 



(3.26) 




(4.1) 



A geostrophic-like model for large Hartmann number flows 



15 



B 




Figure 5. Geometry and magnetic field of a circular duct with fringing magnetic field. 




duct 



axis 



Figure 6. Characteristic surfaces, seen on the (x,y) plane, for half the pipe. The other half 

y < is symmetrical. 



This problem has been initially investigated by Holroyd and Walker [HW78] and later 
by Hua and Walker [HW89]. In the next section 4.1 a two-dimensional numerical calcu- 
lation of the flow is performed using the two-dimensional model described in section 3.1. 
This is somewhat similar to the work done by Hua and Walker [HW89] , but higher values 
of the Hartmann number are investigated here. In section 4.2, an asymptotic model used 
by Holroyd and Walker is re-expressed in terms of the streamfunctions tj) and h instead 
of pressure and electrical potential, and is then solved with a higher accuracy (in effect 
with more terms in a series of Chebyshev polynomials). In this model, it is assumed that 
the non-uniform region is short enough so that the flow follows the characteristic surfaces 
exactly on its length-scale. It will be shown that the numerical two-dimensional model 
and this asymptotic model are in perfect agreement. They both lead to layers of thickness 
Ha" 1 / 4 and it will be shown that the two-dimensional numerical results approach the 
asymptotic results when iJa -1 / 4 is small compared to one. 



4.1. Numerical two-dimensional results 

The two-dimensional equations (3.14) and (3.15) are solved in this section, except for the 
two-dimensional viscous term. This is the last term is equation (3.14) and it has been 
shown in section 3.1 that the boundary or free layers arising from it are always as thin 
as parallel layers, Ha" 1 / 2 , hence cannot be modelled correctly by the two-dimensional 
model. Consequently, this term is removed from the equation and the condition of no-slip 
on solid boundaries is dropped accordingly. These equations are re-written here, in the 
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Figure 7. Typical mesh, with fewer vertices than used in calculations for clarity, for Ha = 10 6 . 
Boundary conditions for velocity and electric current streamfunctions tp and h are shown. 



form that is solved numerically: 

V. ^^.(W x e z )j e z = Ha V x Vh, (4.2) 

v -(™)*- v (t)* v *- '«» 

The symmetry of the problem with respect to the axis y = is exploited and the domain 
of integration is only half of the initial domain: y > 0, — L < x < L, where L is taken 
large enough so that the flow is established in the x direction at these two ends. The 
development length scales as Ha 1 / 2 and this is taken to be the length of the integration 
domain, L = O.bHa 1 ^ 2 . The results show that the flow is independent of x well before 
the ends. 

A problem with this domain is that its length is large for large Hartmann numbers 
and it becomes difficult to build a mesh with a reasonable number of nodes which 
can cope with the short and large length-scales involved. To overcome that difficulty 
a change in x coordinate is made, which maps the long physical domain onto a fic- 
titious domain of length 2. The change in x-coordinate used here is of the form x = 
2.5 sinh (xf sinh _1 (f/a 1/ ' 2 /5)), where Xf is a fictitious x coordinate (see Fig. 7). This 
type of stretching creates a distortion that is more and more pronounced towards the 
ends. The following boundary conditions arc applied: on y = 0, ip = and ft, = 0; on 
y = 0.5, ■0 = 1 and /;, = 0; on Xf = ±1, dijj/dxj = and dh/dx = 0. 

Two types of two-dimensional plots will be shown in this section, none of them involving 
the artificial x f coordinate: one type is a general view of a large part of the domain where 
the xHaT 1 ! 2 will be used while the other type is a limited view of the physical domain, 
between x = —2 and x = 2, where coordinate x will be used. 

The equations are solved using FreeFem+ [PHB+], a free-licence software developed 
at INRIA. This is a two-dimensional finite element solver, coupled to an anisotropic 
unstructured mesh generator that can automatically refine the mesh where solutions 
have large gradients. Triangular, first order elements are used. Typically, after a couple 
of mesh adaptations, the mesh looks like the one displayed on Fig. 7, although the meshes 
used for the calculations contain many more vertices (up to 50, 000). 

It can be seen from Fig. 8 that the large scale structure of the flow and electric cur- 
rent circulation becomes independent of the value of the Hartmann number, at large 
Hartmann number, i.e. more than 10 or so. Between Ha = 10 4 and Ha = 10 6 , some 
differences in position of the iso- , lines may be observed and also of the iso-h lines. It 
is also confirmed from these results that the development length of the flow upstream 
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Figure 8. FreeFem calculations: velocity and electric current streamlines, iso-tf) and iso-/i on 
the left and right respectively for 4 values of the Hartmann number, Ha = 10 2 , 10 4 , 10 6 , 10 8 . 
The vertical coordinate is the y coordinate from y = (centerline of the duct) to y = 0.5, the 
horizontal coordinate is Ka~ x ^ 2 x from Ha~ x ^x = —0.25 to Ha~ 1 ^ 2 x = 0.25. 



and downstream, i.e. along the characteristic surfaces is of order Ha 1 / 2 . However, the 
plots in Fig. 8 do not provide a clear information on the structure of the flow around the 
region of varying magnetic field. This is due to the fact that the axial coordinate is scaled 
with a factor Ha 1 / 2 in order to give a global view of the flow, hence the central region, 
—2 < x < 2 is squeezed enormously in the x direction. The following Fig. 9 shows the 
same results as Fig. 8 when only the central part of the duct is singled out. It can also be 
seen that the flow becomes independent of the Hartmann number, when it is larger than 
10 6 . It can also be seen that, on this length-scale of order unity, the velocity and electric 
current streamlines become closer and closer to the characteristic surfaces (compare with 
figure 4) at large Hartmann number. 

It is also interesting to look at the velocity profiles of the x-component of the velocity 
field across the duct. On Fig. 10, the velocity profile is shown between y = and y = 0.5 
at the axial positions x = — 2 and x = 2 respectively. 

One can see that the velocity profile presents a singularity at infinite Hartmann num- 
ber. This singularity is situated on the characteristic surface which splits in the region 
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Figure 9. FreeFem calculations: velocity and electric current streamlines, iso-^i and iso-h on 
the left and right respectively for 4 values of the Hartmann number, Ha = 10 2 , 10 4 , 10 6 , 10 8 . 
The vertical coordinate is the y coordinate from y = (centerline of the duct) to y = 0.5, the 
horizontal coordinate is x from x — —2 to x = 2. 



of non-uniform magnetic field. Hence, its position is at y = when the velocity profile is 
shown at x = 2, and at y c = 0.5 cos(7r/6) at x = —2. This critical position y c corresponds 
to the y position where the depth of the cylinder in the direction of the magnetic field 
is half the diameter of the duct. On Fig. 11, the velocity profiles are plotted again, in 
logarithmic coordinates, to put in evidence the singularity arising at large Hartmann 
number. On the left (x = —2), only the region y c < y < 0.5 is shown and the horizon- 
tal coordinate used is y — y c . On the right (x = 2), the singularity being situated at 
y = 0, no change of coordinate is made. It can be seen that, at x = 2, the velocity pro- 
file behaves asymptotically like y~ 1 ^ 2 when Ha increases. On x = —2, it looks plausible 
that the velocity behaves as (y — y c )~ 3 ^ 4 when Ha increases: it is not as clear as the 
— 1/2 exponent at x = 2, probably because, around y = y c , the depth varies much more 
rapidly with y than around y = 0. However, there must be a strong relationship between 
these exponents on each side of the change of magnetic field if one believes that, on the 
short length-scale of the change in magnetic field, the streamlines must coincide with the 




Figure 10. FreeFem calculations: velocity profile (x component) at the position x = —2 (left) 
and x = 2 (right), between y = (centerline of the duct) and y = 0.5 (edge of the duct), for 
Ha = \Q 2 , 10 3 , 10 4 , 10 5 , 10 6 , 10 7 , 10 8 , 10 9 , 10 10 . 

characteristic surfaces 

i\x = -2, m (K)] = il>[x = 2,y r (K)\, (4.4) 

where yi(K) is the y position on the left-hand side of the change of magnetic field for a 
given value K of the characteristic function B z /d and where y r {K) is the corresponding 
y position on the right (x = 2). Differentiating expression (4.4) with respect to K , one 
gets: 

— — [x = -2, Vl (K)] = ^^[x = 2,y r (K)]. (4.5) 

The derivatives of yi and y r with respect to K, near the singular characteristic value 
K = 2/1 = 2, are the inverse of the derivatives of K with respect to y, near y = y c and 
y = respectively. From K = B z /d, with d = 0.5(0.25 — y 2 ) 1 ^ 2 , the result: 

1 - 2 

A ~ 0.5(0.25-^)1/2 - o.5(0.25 -y*) 1 / 2 (4 ' 6) 

is obtained, which when differentiated leads to: 

^ = 0.5 yi (0.25 - y 2 )- 3 / 2 , = y r (0.25 - y 2 )^ 2 . (4.7) 

dyi dy r 

This means that, in the neighbourhood of yi = y c and y r = respectively: 

^=8V3, d ^=8y r , (4.8) 

where for the case of dK/dyi the only interesting matter is that it has a finite value which 
is to be considered as uniform near the singularity. In contrast, for the case of dK/dy r a 
first-order expansion is considered since its value at y r = is zero and cannot be inverted 
to feed into equation (4.5). Besides, using (4.6), it can be shown that near the singularity, 
the following relationship holds between yi and y r : 

y?~2V3(yi-y c ). (4.9) 

With (4.8) and (4.9), if it supposed that the the velocity profile in the x direction, dip/dy r , 
on the right-hand side behaves near y r = as u x — ay h r , then equation (4.5) indicates 
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how the velocity should behave on the left-hand side: 

dtp 



dyi 



(x = -2, yi) = V3a 2v / 3(y ( - y c ) 



(4.10) 



This relationship indicate which power law should be found on the left of the singularity 
given that an exponent b is observed on the right. As the exponent b = —1/2 is clearly 
visible on Fig. 11, one can infer from (4.10) that the singularity in the velocity profile on 
the left is of the form (y; — y c )" 3 / 4 . 

Having established the nature of the singularity developing along this particular char- 
acteristic surface which splits as a result of the change of magnetic field, one can also 
examine how this singularity is smoothed out at high but finite Hartmann number. As 
predicted in section 3.1 there is in fact a free shear layer of thickness Ha" 1 ! 4 . This is 
because the change in magnetic field occurs on a typical length scale unity, so when the 
velocity profile is examined at a distance of order unity from the change in magnetic 
field, a layer of thickness Ha -1 / 4 is expected. This can be seen on Fig. 10 and 11 as the 
singularity is smoothed on a shorter and shorter length scale, <5, as the Hartmann num- 
ber is increased. One can actually extract quantitatively the 'smoothing' thickness out 
of these velocity profiles. If we assume that the velocity u x is at its maximum u max on 
a length-scale 5, then the following integrals of the velocity profile can be approximated 



as: 



0.5 



u 4 dy ~ 8u Ti 



0.5 



(4.11) 



io Jo 
as the contribution around the maximum velocity dominates the integrals. It then follows 
from (4.11) that the following value for S is an objective measure of the thickness of the 
layer which can be computed from the profiles: 



5 = 



Jo' 5 u i d y 



1/4 



Jo°' 5 U x d V 



1/8 



(4.12) 



This value of 8 is plotted on Fig. 12 (left). It can be seen that at large Hartmann number, 
5 behaves as Ha" 1 / 4 . This is a practical application of the general analysis in section 3.1. 
indicating that layers of thickness Ha" 1 / 4 can develop along characteristic surfaces when 
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Figure 12. Expression (4.12) is plotted for Hartmann number from 10 2 to 10 12 , based on the 
velocity profile at x = —2 (left) and x = 2 (right). The dashed lines are straight lines of slope 
— 1/4 and —1/6 for comparison. 



a change occurs on a length-scale unity. Here the magnetic field changes on a length 
scale unity and the layer between the stagnant and moving fluids is indeed of thickness 
Ha -1 / 4 . These findings bring an answer to the discussion started in the conclusion of 
the paper by Holroyd and Walker [HW78] . In their conclusion, about the ' nature of the 
velocity profile near the moving fluid/stagnant fluid boundary in the non-uniform region', 
the two authors express different views: 1 J.S.W believes the velocity gradients to be large 
yet Oil) while R.J.H. believes that a shear layer of thickness Oi^Ha -1 / 2 ) might separate 
the moving and stagnant fluid . The answer proposed here is intermediate, with a layer 
of thickness Ha^ 1 ^ 4 . In addition, it is observed that the asymptotic velocity profile itself 
diverges near y c , not only its gradient. 

If one observes now the velocity profile on the right-hand side of the change in magnetic 
field, one can also measure the size of the free-shear layer at y = 0, lying on the dividing 
characteristic surface. Using the same expression (4.12), the thickness of this layer is 
plotted on Fig. 12 (right). It can be seen that, at high Hartmann number, the thickness 
decreases as Ha~ x / 6 . This new exponent arises because the characteristic function K 
shows an cxtremum at y = 0. Indeed, the power law i?a -1 / 4 was derived under the 
assumption that K has a finite gradient G. Here, G is not constant at the scale of the 
layer, but is a linear function of y, G = cz, where c is a constant. When substituting this 
expression for G in equation (3.20), it is possible to carry on a similar scaling analysis as 
in 3.1 to find out that a shear layer of thickness Ha^ 1 / 6 should develop near y = for 
the axial position x = 2. 

4.2. The model of Holroyd and Walker revisited 

The method of eigenfunction expansion used by Holroyd and Walker [HW78] (directly 
adapted itself from a paper by Walker and Ludford [WL74] ) has been followed with two 
changes. First, the problem is solved here using the streamfunctions ip and h, rather the 
the pressure and electric potential p and 0, and secondly, as many as 1000 eigenvectors 
have been used while Holroyd and Walker could only compute 60 at the time of their 
publication. This increase in computer power is necessary to observe the results found in 
the previous section with FreeFem. In this method, the magnetic field is modelled as a 
pure step function and the velocity profiles are calculated at a distance x = 2 upstream 
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Figure 13. Eigenfunction expansions: velocity profile (x component) at the position x = —2 
(left) and x = 2 (right), between y = (centerline of the duct) and y — 0.5 (edge of the duct), 
for.ffa = 10 2 , 1(T, 10*, 10 5 , 10 6 , 10 7 , 10 s , 10 9 , 10 10 . 




Figure 14. Eigenfunction expansions: velocity profile (a; component) at the position x = —2 
(left), between y — y c and y = 0.5 and x — 2 (right), between y = and y — 0.5, for Hartmann 
numbers Ha = 10 6 , 10 7 , 10 s , 10 9 , 10 10 , 10 11 , 10 12 . On the left, the abscissa is the y coordinate 
measured from the singularity y — y c . 



and downstream of the discontinuity, and plotted in linear (Fig. 13) and logarithmic 
(Fig. 14) coordinates. 

The results of this modelling are very similar to those obtained by finite element 
analysis (Fig. 10 and 11). This can be seen as a confirmation that the approach of 
Holroyd and Walker was correct. The results arc not identical though for two reasons. 
First, at Ha = 10 12 and to a lesser extend at Ha = 10 11 , one can see numerical oscillations 
on the right of Fig. 14, showing that 1000 eigenmodes is not sufficient. Secondly, the peak 
values in shear layers are not identical because x = —2 or x = 2 does not have the same 
significance in both models. They are the ends of the gradient region of magnetic field 
in the finite element modelling while they arc exact positions with respect to a sharp 
change in magnetic field intensity in the method of matched expansions. 
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5. The entry problem in rectangular ducts 

The so-called entry problem is another typically relevant issue regarding the two- 
dimensional MHD flow structure. This phenomenon concerns the changes in a duct flow 
when the intensity of an imposed transverse magnetic field varies along the direction of 
the flow, for instance when a duct flow enters (or leaves) a magnet, hence the term 'entry 
problem'. In electrically insulating rectangular ducts, it is observed that the flow remains 
quasi two-dimensional but concentrates near the two parallel walls to the magnetic field 
while becoming very weak in the middle of the duct: because of the no-slip condition at 
the wall, this has been described as an M-shape velocity profile when plotted from one 
parallel wall to the other. An entirely equivalent problem is that of a flow in a uniform 
magnetic field but with a variable cross-section area along the duct. Indeed, only changes 
in the value of the characteristic function B z /d are relevant. Here, we choose to consider 
the case of a uniform depth d — 1 and a variable magnetic field intensity B z as a function 
of x, the longitudinal coordinate along the duct axis. As before z is the coordinate in 
the direction of the magnetic field and y is the other coordinate in the plane of the 
cross-section. 

This problem will be treated using equations (3.14) and (3.15), and our choice leads 
to simplified governing equations: 

Ha [V. (B g Vip)] e z = Ha 2 VB z x V/i + [(V 2 ) 2 ^] e 2 , (5.1) 
[V 2 h] e z = VB Z x Vip. (5.2) 

Then, following our analysis in section cr, the last term of equation (5.1) - of viscous 
origin - will be neglected as it leads to thin parallel layers, too thin for the condition 
of two-dimensional flow to apply. Assuming that its role will be simply to ensure the 
no-slip condition, this term is removed while the no-slip condition is dropped. The set of 
equations becomes: 

[V. (B 2 VV>)] e z = Ha\7B z x V/i, (5.3) 
[V 2 h] e z = VB Z x W- (5.4) 

These equations are solved in section 5.2, for the case of a duct flow with a ramp of 
transverse magnetic field. Before, in section 5.1, a local analytical solution for equations 
(5.3) and (5.4) is derived corresponding to boundary layers analogous to Stommel layers 
for rotating flows. Finally, in section 5.3, an attempt is made to calculate the flow in the 
duct cross-section in a region of axially varying B z : it is compared to and gives support 
to the two-dimensional results. 



5.1. Analytical local solution 

One considers here a region near a side wall, where the gradient of magnetic field B z 
is denoted G. In the Laplacian terms, the normal derivative along y is assumed to be 
dominant compared to the axial x derivative. Equations (5.3) and (5.4) can be written: 
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Since B z and G are functions of x only, these equations are ordinary linear differential 
equations in y. They admit elementary solutions of the form: 



exp I y u Hr) - cxp \- y y w I > (5J) 

explaining that the flow is confined near the lateral walls in a boundary layer of typical 
thickness Ha~ 1 / 2 G~ 1 , when the scale for magnetic field strength is chosen locally {B z = 
1). As long as the gradient of magnetic field G is very small compared to one, the two- 
dimensional structure of the layer is guaranteed and its velocity and electric current 
density must be exponential. 

5.2. Numerical two-dimensional results 

A duct of uniform dimensionless depth unity, width equal to 2, and length 8 is considered. 
The upstream and downstream region occupy each one fourth of the length of the duct 
and have a uniform transverse magnetic field, of intensity 1 — 2G and 1 + 2G respectively, 
while the middle part is submitted to a uniform gradient G, B z = 1 + Gx. On the lateral 
walls, y = ±1, the electrically insulating condition can be written h = and the condition 
of impermeable walls is written i\) = on the lower wall at y = — 1 and i\) = 1 at y = If. 
At the ends of our duct, we apply the condition that h and ip have no normal gradient. 
This is consistent with the fact the flow is nearly fully established at the ends. 

The streamlines of the computed solution are shown on Fig. 15 for Ha = 10 5 and three 
values of the gradient of magnetic field applied, G = 0.05, G = 0.1 and G = 0.2. It can 
be seen that the side layers become thicker as G decreases. In the real three-dimensional 
case, it is expected that a parallel layer will exist, as a sublayer within this thicker side 
layer. The exponential velocity distribution (5.7) has been compared with the calculated 
distribution with FreeFem: the agreement is excellent. 

5.3. Local solution in a cross-section 

As we have seen in section 5.1, the solution of the duct flow with a gradient of magnetic 
field is local: it is sufficient to know the local value of the magnetic field B z and the local 
value of the gradient G to determine completely the boundary layer solution, providing B z 
and G do not vary significantly on an axial length-scale comparable to the layer thickness. 
As a consequence, it is expected that the three-dimensional flow can be computed locally 
in a cross section, within some simplifying assumptions. 

Let us assume that the flow will be mainly streamwise, and that its streamwise gradient 
is small compared to its gradient in the cross-section. Navier-stokes equation can be 
written in the streamwise direction: 

-^- + Ha 2 j y B z + V 2 s u x = 0, (5.8) 

where V| is the Laplacian operator in the cross-section. This equation (5.8) will be one 
of our three equations within the cross-section, the u x equation, and we shall now derive 
equations for <p and dp/dx. The unknown j y in (5.8) can be expressed as a function of <f> 
and u x according to Ohm's law: j y — —d(j)/dy — u x B z . 

The equation for (f> is obtained as usual by taking the divergence of Ohm's law: 

V 2 = B.Vxu, (5.9) 

f The condition -0 = 1 can be changed for an arbitrary constant without affecting the nature 
of the solution as we are treating a linear problem. 
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Figure 15. Velocity streamlines (left) and electric current streamlines (right) resulting from the 
two-dimensional numerical analysis for Ha = 10 5 and three values of G, G = 0.05 on the top, 
G — 0.1 in the middle and G = 0.2 at the bottom. The vertical coordinate is the y coordinate 
from y = (centerline of the duct) to y = 1, the horizontal coordinate is x from x — —4 to to 
x = 4. The symmetrical part of the cavity (— 1 < y < 0) is not represented here. 

With a dominant velocity cc-component, and with gradients predominantly in the cross- 
section, this equation can be written: 

V^ = -B,^, (5.10) 
Oy 

An equation for dp/dx is a less standard object. It is obtained through a few steps. 
First, we consider the curl of the Navicr-Stokcs equation: 

Ha 2 B z ^- - Ha 2 Gj x e z + V 2 (V x u) = 0. (5.11) 

oz 

As a next step, we shall write the x-component of the curl of equation (5.11): 

Ha 2 Bz d(yxih _ Ha 2 G d _^L _ (V ^ U;c = 0, (5.12) 
oz oy 

To make progress, we need to derive the curl of Ohm's law: 

n 

V xj = B z ^-u x Ge z , (5.13) 
oz 

where e z is the unit vector in the z direction. Substituting into (5.12), and assuming 
again that gradients are essentially due to variations with the cross-section, we get: 



Ha 2 Bj^ - Ha 2 G 2 u x - (V 2 s j 2 u x = 0, (5.14) 



On the other hand, one can derive a similar equation by taking directly the cross-section 
Laplacian of equation (5.8): 

-V|^ + Ha 2 B 2 W 2 Jy - {V 2 s ) 2 u x = 0. (5.15) 



26 



Thierry Alboussiere 



Taking the curl of equation (5.13) and considering its y component leads to: 

~V 2 sJy =B^, (5.16) 

as variations along x are neglected compared to variations within the cross-section. Sub- 
stituting into (5.15) leads to: 

-V|g-ffa 2 B^-(V|)V = 0. (5.17) 
The difference between equation (5.14) and (5.17) gives an equation for dp/ dx: 

-V|^ = Ha 2 G 2 u x . (5.18) 
ox 

In summary, we have a system of three equations for three unknowns: 

V 2 s u x = ^- + Ha 2 B z ^ + Ha 2 u x Bl (5.19) 
ox ay 

V|0 = -B^, (5.20) 
dy 

V| g = -Ha 2 G 2 u x . (5.21) 

This set of equations is solved for the same geometry as considered in section (5.2). The 
cross-section has dimension 1 along the magnetic field direction and 2 in the transverse 
direction. Because of the symmetries of the configuration, we solve this problem in a 
quarter of the cross-section. 

Finally, we need to specify boundary conditions for u x , <j) and dp/dx. The conditions 
on u x and </> are obvious, with u x vanishing on the walls, and dcf>/dn — as electrically 
insulating boundaries are considered. Only approximate boundary conditions can be 
found for dp/dx; the assumption that this variable is going to be nearly independent of 
z is made. This assumption is true within the core, it is also quite true within Hartmann 
layers as they are so thin that pressure cannot change significantly across them, and we 
also assume that three-dimensional parallel layers (-Ha -1 / 2 ) are thin enough so that they 
cannot cause large z-dependence of dp/dx. Then, integrating equation (5.21) across the 
whole cross-section tells us that, along parallel walls, the normal derivative of dp/dx is 
proportional to the total volume flow rate though the duct, Q: 

while this derivative vanishes on Hartmann layer walls. An arbitrary value Q = 2 will be 
used in the calculations so that the mean velocity is unity. 

The distribution of axial velocity, electric potential and axial pressure gradient was 
calculated with FrccFcmH — h, an improved version of FreeFem+ where it is possible to 
used triangular second order finite elements. This proved necessary as Hartmann layers 
were resolved for Hartmann number up to 10 . In this cross-section modelling, Hartmann 
layers have to be very well resolved as the electric current flowing through them is a key 
ingredient in the establishment of Ha~ x l 2 G~ x layers. 

Some results are shown in Fig. 17, for a fixed gradient G = 0.2, and increasing Hart- 
mann numbers, from 10 to 10 5 . As expected, shear layers of thickness Ha~ 1 / 2 G~ 1 develop 
and parallel layers of thickness Ha~ x l 2 arc sublayers: their role is indeed to bring the 
velocity to zero at the wall. The velocity profile at z = is compared with the analytical 




Figure 16. Numerical calculation of the flow in the cross-section of a rectangular duct. Using 
symmetries, only a quarter of the cross-section is analyzed. Isovalues of the streamwise velocity 
(left) for Ha = 10, 10 2 , 10 3 , 10 4 and 10 5 , with G = 0.2. Velocity profile at z = (right) for the 
same parameters, compared to the analytical solution (5.23). 
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exponential solution. By symmetry, using expressions (5.7) and for a total volume flow 
rate equal to 2 (Q = 2), the analytical expression (dashed line in Fig. 17) is: 



The agreement is very good, the small departure being due to the flow deficit in parallel 
layers at large Hartmann number and to the deficit in Hartmann layers at smaller Hart- 
mann numbers. This shows that our local problem in a cross-section (equations (5.19), 
(5.20) and (5.21)) is a good representation of MHD flows for entry problems. It is also a 
direct confirmation that the results of two-dimensional modelling can be confirmed from 
a model where Hartmann layers are effectively calculated. 

6. Concluding remarks 

The structure of MHD flows under large Hartmann numbers has been investigated 
for an arbitrary geometry of the container and magnetic field distribution. A set of two 
equations governing the two-dimensional mass flux streamfunction ip and electric charge 
flux streamfunction h has been used, however it would be equivalent to consider other 
variables such as pressure p and electric potential cf> as in [HW78] . The reason why mass 
and electric fluxes are preferred is that this approach shows clearly which properties of 
Hartmann and Ekman layer play a key role in the two-dimensional equations (see section 
3). In the case of rotating flows, the key property of Ekman layers is to carry a mass flux 
in the direction perpendicular to the core velocity. In the case of MHD flows, the key 
property of Hartmann layers is to carry an electric current in the direction perpendicular 
to the core velocity. 

The originality of the present work consists largely in the scaling analysis of these MHD 
two-dimensional equations. It has been shown that shear layers of thickness Ha -1 / 4 can 
develop. In addition, the analogy between rotating flows and MHD flows has been further 
developed, so that for instance it is now clear that the analogue of Stommel layers are 
those analyzed by Walker and Ludford [WL72] of thickness of order Ha~ 1 / 2 G~ 1 . 

It may be that some of the results presented in this paper, especially those concerning 
layers of thickness Ha^ 1 ^, will never realistically apply. In industrial MHD situations, 
one may envisage Hartmann numbers up to 10 5 , with 10 Tesla and 10 cm length-scale for 
a liquid metal. However, even in this case, it may be interesting to know what the flow 
structure would be at higher Hartmann numbers, as a tendency towards this structure will 
already exist at lower Hartmann numbers. In the astrophysical and geophysical context, 
it is easy to find higher estimates of Hartmann numbers. For instance, in the liquid 
Earth core, one can reach Ha ~ 10 7 or Ha ~ 10 8 depending on whether this is based 
on a poloidal (B ~ 5 10~ 4 ) or a toroidal (B ~ 5 10~ 3 ) estimate for the magnetic field 
intensity (see [CBNM02]). Nevertheless, rotating effects, inertial effects and transport of 
magnetic field are so important that the Hartmann number is not the key dimensionlcss 
parameter. 

A number of possible extensions of this work can be envisaged. First, in this paper, the 
use of the straight magnetic lines approximation has been made and a more accurate two- 
dimensional model could be obtained by taking into account of the true magnetic field 
distribution. Although the results would be changed quantitatively, it is not expected 
that the flow structure would be altered significantly. Secondly, the present analysis can 
be extended to analyse free-surface flows. The position of the free surface then becomes 




(5.23) 
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an additional unknown and the nature of Hartmann layers changes : they become much 
less active in the sense that the electric current flowing along them is drastically reduced 
compared to a wall-bounded Hartmann layer. Shear layers of a new type are expected 
as the balance in our equation (3.18) will be changed. A third possible extension would 
be to consider the combined effect of rotation and magnetic field. This is typical of 
astrophysical or planetary dynamics. When rotation and magnetic field have an identical 
direction, there is certainly a predominantly two-dimensional flow. In the general case of 
different orientations, it is far from obvious that such a two-dimensional dynamics will 
develop, unless one effect (rotation or magnetic field) is dominant compared to the other 
one. 

The analogy with rotating flows has been used as an example to follow in the MHD 
case. However, if one takes the opposite view and wonder what should be the analogous 
of Ha" 1 / 4 layers, one can go back to equation (3.26) for the homogeneous model. If one 
considers a change on a length-scale unity along a gcostrophic contour, one may wonder 
what the thickness 6 of the shear layer along that contour should be. In Oceans, as [3 
is small, the answer is that it should be of order iS 1 / 4 /^ 1 / 2 and that it is the result 
of the competition between topographic effects and Ekman pumping effects. This layer 
is thicker than the E 1 ^ Stewartson layer, which corresponds to the balance between 
bulk viscous effects and Ekman pumping. Are there shear layers, lying along geostrophic 
contours, of thickness of order E 11 / 4 ^ -1 / 2 ? 
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